##安装"DESeq2"##
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("DESeq2")

##数据处理##
library("DESeq2")
df <- read.csv("E:\\class7\\mRNA_exprSet.csv")
rownames(df)=df[,1]   #将第一列作为行名
df <- df[,-1]            #将第一列删除
df1 <- df[rowMeans(df)>1,]  #去除表达量过低的基因
group <- substr(colnames(df1),14,15)  
table(group)            #分组信息
colData <- data.frame(row.names=colnames(df1), group)  #样本信息

##差异表达分析##
dds <- DESeqDataSetFromMatrix(countData = df1, colData = colData, design = ~ group)
head(dds)

##标准化##
dds1 <- DESeq(dds) 
res <- results(dds1)
summary(res)   #查看结果

##筛选表达量显著上调/下调的基因##
# res格式转化：用data.frame转化为表格形式
res1 <- data.frame(res, stringsAsFactors = FALSE, check.names = FALSE)
# 依次按照pvalue值log2FoldChange值进行排序
res1 <- res1[order(res1$pvalue, res1$log2FoldChange, decreasing = c(FALSE, TRUE)), ]
# 获取padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因
res1_up<- res1[which(res1$log2FoldChange >= 1 & res1$pvalue < 0.05),]      # 表达量显著上升的基因
res1_down<- res1[which(res1$log2FoldChange <= -1 & res1$pvalue < 0.05),]    # 表达量显著下降的基因
res1_total <- rbind(res1_up,res1_down)

##可视化：绘制火山图##
library("ggplot2")
library("ggrepel")
genes<- res1
genes$color <- ifelse(genes$padj<0.05 & abs(genes$log2FoldChange)>= 1,
                      ifelse(genes$log2FoldChange > 1,'red','blue'),'gray')  #红色表示显著上调基因，蓝色表示显著下调基因，其余为灰色
color <- c(red = "red",gray = "gray",blue = "blue")

p <- ggplot(
  # 指定数据、映射、颜色
  genes, aes(log2FoldChange, -log10(padj), col = color)) +  
  geom_point() +
  theme_bw() +
  scale_color_manual(values = color) +
  # 辅助线
  labs(x="log2 (fold change)",y="-log10 (q-value)") +
  geom_hline(yintercept = -log10(0.05), lty=4,col="grey",lwd=0.6) +
  geom_vline(xintercept = c(-1, 1), lty=4,col="grey",lwd=0.6) +
  # 图例
  theme(legend.position = "none",
        panel.grid=element_blank(),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 14)+
          # 注释
          geom_text_repel(
            data = subset(genes, padj < 1e-100 & abs(genes$log2FoldChange) >= 10),
            aes(label = rownames(genes)),
            size = 5,
            box.padding = unit(0.35, "lines"),
            point.padding = unit(0.3, "lines"))
  )
pdf(file="E:\\class7\\volcano.pdf")
p
dev.off()
